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Abstract 

We consider the general character of the spatial distribution of a population that grows through 
reproduction and subsequent local resettlement of new population members. We present several 
simple one and two-dimensional point placement models to illustrate possible generic behavior of 
these distributions. We show, numerically and analytically, that these models all lead to multifrac- 
tal spatial distributions of population. Additionally, we make qualitative links between our models 
and the example of the Earth at Night image 1], showing the Earth's nighttime man-made light 
as seen from space. The Earth at Night data suffer from saturation of the sensing photodetectors 
at high brightness ('clipping'), and we account for how this influences the determined dimension 
spectrum of the light intensity distribution. 
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FIG. 1: A version of the Earth at Night image []|. This is a composite image created with data 
from the Defense Meteorological Satellite Program of the Earth's nighttime man-made lights as 
seen from space. For the analysis, a grayscale 30000x15000 pixel image was used. The image pixels 
have intensity values ranging from to 100. 

I. INTRODUCTION 

Growing populations exist in many different areas of interest. Although the most obvious 
examples are biological in nature (e.g., bacterial cultures, the human population), one can 
also consider the growth of technological development or the growth of urban infrastructure 
as more abstract examples. While the specific details of these systems can be quite different, 
the examples stated here all share the characteristic that their populations grow on some 
background space. A reasonable question to ask then is what types of spatial distributions 
result in these growing systems. Figure ^ is a version of the popular Earth at Night image 
(EaN) jl|. This is a composite image taken by orbiting satellites that shows the Earth's 
nighttime man-made lights as seen from space. The light intensities are brightest in areas 
that are known to be highly populated and developed, and thus the image can be crudely 
thought of as representing the spatial distribution of some combined measure of technological 
development and human population density. We see a very heterogeneous distribution, with 
areas of very high intensities as well of almost no light at all. 

A useful characterization of the heterogeneity of a distribution is the distribution's di- 
mension spectrum Dg. A grid of hypercubes, or boxes, of a fixed size with edge length e is 
used to cover the distribution, where the dimension of the hypercubes corresponds to the 
dimension of the underlying space (e.g., squares for two dimensions and line segments for 
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one dimension). Dq is then calculated via 

" ln(l/e) ' ^ ^ 

where 

N{e) 

me) = Y.f^l- (2) 

i=l 

Here the sum is taken over A^(e) non-empty grid boxes and /ij is the fraction of the partic- 
ular quantity of interest (assumed to be nonnegative, e.g., light intensity in Fig. [TJ that is 
contained in box i. 

The parameter g (g > in this paper) can be varied continuously to vary the influence 
of high and low /Xj boxes, with larger values of q emphasizing the boxes with larger /ij's. In 
particular, when g = 1 we have the information dimension Di which, applying L'Hospital's 
Rule to (H)) and (0), is given by 

jj ^.^l^^^lJ^^ (3) 

Ine ^ ^ 

A distribution is fractal if it possesses a fractional (non-integer) dimension (i.e., if Dq is not 
an integer for some q). If Dq depends on q, then the distribution is said to be multifractal. 

Generally, when measuring the dimension of a distribution of numerical and/or experi- 
mental data, there is a maximum resolution observable. In the EaN image, for example, the 
image's individual pixels are the smallest boxes that can be used to cover the distribution in 
the calculation of Dq. In such cases the e — > limit cannot be considered; instead one looks 
for a scaling range of e where the quantity (1 — q)^^ ln/(g, e) varies approximately linearly 
with ln(l/e), from which a dimension can be extracted as the slope of a straight line fitted 
to data in the e scaling range. We illustrate this process by calculating the information 
dimension Di of North America (6000x3000 pixels) in the EaN image. Figure El shows a 
plot of Y!i=i A^j In /^i versus In e (with e in units of pixels) for the image, where each /ij is 
the fraction of the total intensity in box i. The good evidence of linearity in this plot over 
a scaling range in e of order ~ 50 indicates that a fractal description makes sense. The 
solid line in Fig. El has a slope of 1.65 ± 0.02, which we take as the value of Di (the error 
indicates the uncertainty involved in extracting the slope of the scaling region). Thus the 
light distribution over North America in the EaN is fractal. (Other regions of the EaN also 
exhibit fractality, with different values of Di.) In Sec. Illll we demonstrate that the EaN light 
distribution is also multifractal. 
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FIG. 2: Plot for calculating the information dimension Di [see Eq. Q] for North America in the 
Earth at Night image. The quantity J2i^=i^ f^'i l^Mi is plotted versus Ine (open circles), where e is in 
units of pixels, resulting in a linear scaling region, the slope of which is a finite-scale approximation 
of Di. The solid line has a slope of 1.65 it 0.02, where the error value indicates the uncertainty 
involved in extracting the slope of the scaling region. 

Yook et al. ^ have recently reported that the world's human population, as well as the 
population of internet routers and autonomous systems, are fractally distributed. Addition- 
ally, several studies have demonstrated various fractal properties of urban 
settlements. We were thus led to consider the possibility of common, underlying mechanisms 
shared by growing populations that lead to fractal distributions. 

One characteristic found in many growing populations is the existence of some general- 
ized reproduction process, where existing members of the population generate more mem- 
bers. In addition to biological reproduction, in some non-biological systems, concentrated 
populations can encourage the creation of more members. For example, areas that are tech- 
nologically highly developed are likely to stimulate more development than areas that are 
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not as technologically developed. 

Generally in such systems, new members of a population become situated near an existing 
member, but not at the exact same location. We model this process as a local "resettlement" , 
whereby an existing member of the population produces one or more new members, who 
then resettle to a nearby location. 

In this paper we demonstrate that these two ingredients, reproduction and local reset- 
tlement, can lead to multifractal spatial distributions for growing populations. We present 
several point placement models in one and two dimensions which implement these ideas and 
show, both numerically and analytically, that they lead to multifractal distributions. 

One might suspect that fractal population densities arise due to strong inhomogeneities 
of the underlying space (e.g., inhomogeneities in the distribution of land fertility, natural 
resources, etc.). Our models initially (Sec. |n} employ spaces that have no inhomogeneity, 
thus demonstrating that fractal population distributions can occur without inhomogeneity 
of the background space on which the population grows. We also consider generalizations 
to our models that include geographical inhomogeneities (Sec. IIV|1 and demonstrate that, 
for the particular forms of the inhomogeneities that we employ, the multifractality of the 
distribution is unaffected. 

In Sec. |n] we present our point placement models along with numerical results. (We 
emphasize that the set of models we introduce in Sec. |n]is not exhaustive, and many other 
similar models could be conceived.) Section IIIII considers situations where a multifractal 
distribution is sensed by an instrument that saturates at a maximum measurable value 
(clipping). We show that such clipping occurs in the EaN data, affecting the determined 
Dg. In Sec. IIVI we investigate the effect of adding geographical inhomogeneity to the under- 
lying space of our models. In Sec. |3 we present analytical results for some of our models. 
We summarize our findings and conclude in Sec. I VII Appendix ^provides details and back- 
ground on the analytical calculation of the spectra of fractal dimensions for some of our 
models. Appendix El contains the proof of a theorem used in Appendix El (this theorem 
should be of very general utility, not restricted to our specific models, in analysis of the 
spectra of fractal dimensions). 
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t = 3 



t = 4 



FIG. 3: Three representative steps of Models 1 and 3. 
II. MODELS AND RESULTS 

A. 1-D Random Interval (Model 1) 

We begin with a very simple model that places points on the circumference of a circle. (For 
convenience we consider a circumference of unit length.) Our initial state is a single point on 
the circle's circumference. As the population (number of points on the circumference of the 
circle) grows, we consider the intervals of the circle's circumference between two adjacent 
points. At each subsequent discrete time step (each 'generation') labeled t, a new point is 
'born' in the middle of each pre-existing interval with the same given fixed birth probability 
p, with each new point bisecting its corresponding target interval. Several representative 
steps of this construction are illustrated in Fig. IHl 

In this way, every interval has an equal probability to generate or attract new points 
at the same rate p. This also implies that regions with a high density of points (i.e., with 
a greater number of intervals) are more likely to attract new members than areas that are 
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sparsely populated, creating a "rich get richer" phenomenon. We will see that all the models 
presented in this paper have this general characteristic. The simple rule of placing each new 
point in the center of an interval allows us to obtain analytical results, as will be shown in 
Sec. EH 

For p = 1 each existing interval is bisected by a point at each time step, resulting in 
a uniformly distributed set of points on the circle circumference for all t. Thus, for more 
interesting point distributions, we turn our attention to models with p < 1. As an example, 
we numerically examine the p = 0.5 case. Figure |3] shows semilog plots of point location 
histograms, along the circle circumference, with different bin sizes Ax for one realization 
of our point placement scheme with p = 0.5, at time t = 35. We observe a roughening of 
the plot as smaller histogram bins are used, thus giving a sense of the heterogeneity of the 
distribution. To quantify this heterogeneity, we calculate the dimension spectrum Dq of the 
distribution. We cover the distribution of points with intervals of a fixed length e and assign 
to each interval i a measure fii equal to the fraction of the total points contained in i. Then, 
referring to Eq. (jT}, we plot the quantity (1 — g)~^ln/(g, e) versus ln(l/e) for a range of e, 
so that the slope of the graph in the scaling range gives us Dg. 

For q > 0.5 the plots exhibit a linear regime that allows for a relatively unambiguous 
determination of the dimension. However, for smaller q there tend to be shorter scaling 
regions, which makes the extraction of a reliable Dg more error prone. This is due to the 
fact that the distribution is made up of a finite collection of points. As the grid covering 
the distribution is made finer, grid intervals in the more sparsely populated areas begin to 
contain either or 1 points. Furthermore, if e is made small enough, eventually all the points 
in the distribution are covered by different boxes, and the dimension for the distribution is 
simply the dimension of a collection of individual points, i.e., zero. Thus, in order to get a 
meaningful dimension out of a finite distribution of points, a certain degree of coarseness of 
the covering grid is necessary. For small values of q, Dq relies heavily on boxes containing 
small amounts of the measure. Therefore, as e is decreased, a scale is quickly reached 
where individual points in sparsely populated areas are covered by separate boxes, thereby 
destroying the overall scaling behavior. For larger q, the dense regions have a large enough 
influence on Dg that the effects of the finiteness of the distribution are not observed until 
very small scales are reached, allowing for a relatively larger scaling range. 

Figure El is a plot of the numerically determined average Dq (open circles) for 20 real- 
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FIG. 4: Semilog plots of histograms of point locations on the unit length circumference of a circle 
for a point distribution (~ 10^ points) generated by Model 1 with p = 0.5, at time t = 35, for 
various bin sizes Ax. a) Ax = 1/16, b) Ax = 1/32, c) Ax = 1/64, d) Ax = 1/128. 

izations of point distributions generated by Model 1 with p = 0.5 at time t = 35. (The 
number of points in each realization ranges from 1.2 X 10^ to 2.8 x 10^ points.) The solid 
line represents the analytical result for Dq derived in Sec. IVAl fEq. |T!?|) with p = 0.5). Both 
the numerical and analytical results demonstrate that this model generates multifractal 
distributions. 

Notice that our analytical value for the box-counting dimension Dq is one. Unlike Dg for 
q positive, Dq depends only on the limiting set of points and not their distribution. Our 
result that -Dq = 1 indicates that, although the point distribution is multifractal, the set 
filled by the points in Model 1 as time t goes to infinity is not fractal. In fact, the entire 
circumference is filled with probability one, which can be easily seen as follows. Since each 
interval has a probability (1 —p) of not being bisected at each time step, the probability for 
an interval not to have been bisected by time t is simply (1 — p)*, which approaches zero as t 
approaches infinity. Therefore no empty intervals remain as t ^ oo. One can similarly argue 
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for the other models in this paper (except Model 7) that the limiting sets are not fractal; 
the heterogeneous point distributions they yield are reflected by the fractional values of Dq 
for g > 0. 

We can also consider Model 1 in the limit where p ^ 0. We think of this limit as the 
situation where at most one interval is chosen to be bisected from all available intervals 
at each time step. This allows us to reformulate the point placement scheme as follows. 
We begin with the same initial condition as the original model, which is a single point on 
the circle circumference. Then at each subsequent discrete time step we place a new point 
according to the following prescription: (a) a target point, or "parent", is chosen from all of 
the preexisting points with equal probability; (b) the new point, or "child", is placed in the 
middle of one of the two empty intervals adjacent to the target point with equal probability. 
Steps (a) and (b) are repeated until the desired system size is reached. 

To show that this point placement scheme coincides with the p ^ limit of Model 1, 
we plot the numerically determined average Dg for 20 realizations of 10^ point distributions 
generated by it in Figure El (open triangles) and observe that the data agree with the p — >■ 
limit of the analytical result for Dq (dashed line) calculated in Sec. IV Al [Eq. ()15|) ]. 

Thus we see that Model 1 can create point distributions with Dg ranging from the trivial 
Dq = 1 curve, corresponding to p = 1, to the Dq curve associated with p — > 0. 

B. 2-D Square (Model 2) 

Model 2 is a two dimensional analogue of Model 1, where points are placed on the unit 
square, (0<x<l, 0<y<l), instead of the circumference of a circle. Beginning with a 
single point at (0.5,0.5), we consider reproduction by a process mimicking cell division. At 
each discrete time step, each preexisting point, or "parent", has a probability p of dividing 
into four points, or "children", where these four points disperse to occupy (resettle) the 
centers of the four equal size squares obtained by partitioning the original square associated 
with the parent point; see Fig. IHl 

Note that this construction can be specified without reference to points, but only to 
the squares they occupy: At each discrete time step, each square is subdivided into four 
equal new squares with probability p. After the desired iterations of this area partitioning 
procedure, the corresponding population distribution can be produced by simply placing a 
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FIG. 5: Plots of Dq for Model 1 with p = 0.5 and p ^ 0. The open circles {p = 0.5) and triangles 
{p 0) are the numerical values calculated from the slope of a line fit to (1 — q)^^ lnl{q, e) versus 
ln(l/e) [see Eq. (^], averaged over 20 realizations of point distributions (except for g = 1, in which 
case the line is fit to J2i^=i^ t^i^^t^i versus Ine [see Eq. ©]), with the error bars indicating the 
sample standard deviation for each value (in the p = 0.5 case the error bars are all smaller than 
the open circles and we omit them). The solid line is the analytical result of Eq. (|13j) with p = 0.5, 
while the dashed line corresponds to Eq. (|15j) . As an example the inset shows the determination 
of Di for one point distribution. For p = 0.5, the distributions are calculated at time t = 35 and 
have between 1.2 x 10^ and 2.8 x 10^ points. For p ^ the numerical results were obtained for 
distributions of 10^ points. 
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FIG. 6: Illustration of the four-fold "cell division" step in Model 2. 

point in the center of each square. 

Since p = 1 results in a homogeneous planar distributions of points [Dq = 2), we again 
consider p < 1 for more interesting point distributions. Choosing p = 0.5, we plot in 
Fig. [7| a distribution of points generated by this model at time t = 15. The points form 
a heterogeneous distribution with dense clusters in some regions and relatively sparsely 
populated areas in other regions. Figure |H1 is a plot of the numerically determined average 
Dq (open circles) for 20 realizations of point distributions generated by Model 2, with t = 15 
and p = 0.5. The solid line represents the analytical result for Dq derived in Sec. IV Bl 

As we did for Model 1, we can explore the p — > limit by considering the construction of 
the point distribution one point at a time. The appropriate process is as follows: At each 
time step a "parent" point is chosen out of all preexisting points with equal probability and 
produces four "children" points according to the four-fold cell division process described 
above. The open triangles in Fig. |H1 are the average Dq of 20 independent 10^ point dis- 
tributions generated by this construction, while the dashed line is the p —>■ limit of the 
theoretical result (Sec. IV B|) . Thus, Model 2 creates point distributions with Dq ranging 
between Dq = 2 for p = 1 and the Dq associated with p —>■ 0, and hence, other than the 
p = 1 case, results in multifractal point distributions. 



C. 2-D Triangle (Model 3) 

Model 3 is a slight variation of Model 2. Here triangles, as opposed to squares, are the 
basis of the construction. We begin with a point in the middle of a single equilateral triangle 
with unit edge length. At each time step each preexisting triangle (parent) is divided into 
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FIG. 7: A plot of a point distribution on the unit square generated by Model 2 with t = 15 and 
p = 0.5. 



four identical equilateral triangles (children) with probability p, with the edge lengths of the 
child triangles equal to half the edge length of the original parent triangle (see Fig. ^ and 
with points placed in the centers of the child triangles. We show in Sec. IV Bl that Model 3 
has the identical analytical expression for Dq as Model 2, and thus also results in multifractal 
distributions. 

D. 2-D Random Square (Model 4) 

Model 4 is a two dimensional point placement scheme similar to the p ^ limit construc- 
tion of Model 2, where points are placed on the unit square, (0 < a; < 1, < y < 1). We 
begin with a single point at (0.5, 0.5). Imagining vertical and horizontal lines drawn through 
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FIG. 8: Plot of Dq for Model 2 with p = 0.5 and p ^ 0. The open circles {p = 0.5) and triangles 
(p 0) are the average numerical values of Dq (see caption in Fig. ISJ obtained from 20 realizations 
of point distributions generated by this model and the error bars are the sample standard deviations 
for each value (in the p = 0.5 case the error bars are all comparable in size to the open circles and 
we omit them) . The solid line is the analytical result of Eq. (|21|) and the dashed line is Eq. (|23() . 
For p = 0.5, the distributions are calculated at time t = 15 and have between 1.2 x 10^ and 2.9 x 10^ 
points. For p — > the numerical results were obtained for distributions of 10^ points. 




FIG. 9: Illustration of the "reproductive" step in Model 3. 
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FIG. 10: Three representative steps of Models 4 and 6. 

this point, we see that the point is at the vertex of four neighboring squares. At the next 
time step, we choose one of the squares at random with equal probabihty, and place a point 
in its middle. We then divide this chosen square into four smaller squares by horizontal 
and vertical lines through the newly added point. Then we continue this process: At each 
subsequent discrete time step we pick a target point (parent) from all preexisting points 
with equal probability and add a new point (child) to the center of one of the target point's 
four adjacent squares, also chosen with equal probability. Several representative steps of 
this construction are illustrated in Fig. ^| 

Notice that some squares have two potential parent points at their corners and some 
squares have only one. The former squares are twice as likely to be subdivided as the latter 
squares, unlike the p ^ limit of Model 2, in which all squares are equally likely to be 
subdivided. 

Figure ITT] shows a distribution generated by this model's point placement scheme, where 
we have successively magnified dense regions of the distribution to illustrate the large dif- 
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FIG. 11: A plot of a 4 X 10^ point distribution on the unit square generated by Model 4. Dense 
regions are magnified to illustrate the heterogeneity of the point densities in the distribution at 
different scales. 

ferences in densities that result at various length scales. Figure ^1 is a plot of the average 
Dq of 20 independent distributions generated by this model, where the solid line is the an- 
alytical result (jTTj) derived in Sec. IV (^1 Both the numerical and analytical results confirm 
that Model 4 generates multifractal distributions. 

E. 1-D Larger Interval (Model 5) 

We return now to one dimension with a model similar to the construction of the p — > 
limit of Model 1, but with a different local resettlement rule. Referring to the p — point 
placement prescription of Model 1, the initial state and Step (a) are identical while Step (b) 
is modified. Instead of a new point being randomly placed in one of the adjacent intervals 
of the chosen target point with equal probability, the new point is placed in the larger of 
the two intervals adjacent to the target point, and, if the two intervals are equal in size, one 
is chosen at random. 

In this way, the new points seek out more sparsely populated local intervals to settle 
within. This strategy is reasonable in a context such as human resettlement, where the set- 
tlers may try to combine the convenience of being close to others with the possible advantages 
associated with more space. We expect this resettlement scheme to yield distributions that 
are more homogeneous than those created by the p limit of Model 1, and hence yield 
greater Dq. 

We numerically obtain Dq and find the dimension spectrum to be multifractal, having 
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FIG. 12: Plot of Dq for Model 4. The open circles are the average numerical values of Dq (see 
caption in Fig. [SJ obtained from 20 realizations of 10^ point distributions generated by this model 
and the error bars are the sample standard deviations for each value. The solid line is the analytical 
result of Eq. (|H7|l . 

the same qualitative shape as the p ^ result for Dq for Model 1. However, due to the 
more uniform nature of the distributions created by this model, the specific values of Dq 
are larger at each q (e.g., Di = 0.84 ± 0.01 [averaged over 20 distributions of 10^ points, 
where the error is the sample standard deviation] for this model versus our analytical result 
Di ~ 0.72 for the p limit of Model 1 [Eq. E]). 

F. 2-D Sparse Square (Model 6) 

Model 6 is a two dimensional analogue of Model 5. The initial state and Step (a) of 
its construction are the same as in Model 4. Our aim is to extend Step (b) of Model 5, 
the choosing of the larger interval for settlement, to two dimensions, and, for that purpose, 
we make the following observation. In the square-based two dimensional point placement 
scheme of Model 4, there are two different types of squares which we refer to as type I and 
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FIG. 13: Illustration of type I (labeled I) and type II (labeled II) squares, where black dots indicate 
populated vertices. 

type II squares. A type I square only has one of its vertices occupied by a point, while a 
type II square has two opposite vertices occupied by points. (See Fig. UHlfor an illustration 
of type I and type II squares.) Thus a type I square can be regarded as more sparsely 
populated than a type II square of equal size. With this in mind. Step (b) for Model 6 is 
to place the new node in the middle of the largest square adjacent to the target point. In 
case of a tie, choose a type I square over a type II square, with any subsequent ties being 
resolved by an equal probability random choice from the remaining candidate squares. 

We plot a distribution generated by this model in Fig. We find this distribution to 
be multifractal, with a Dg that is similar in shape to, but everywhere larger than, the Dq 
of Model 4 (except for q = where the two coincide with Dq = 2). In particular, the 
information dimension is Di = 1.30 ±0.02 (averaged over 20 distributions of 10^ points), as 
compared to our analytical result of Di ^ 1.05 for Model 4 (Eq. IHH]) . 
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FIG. 14: A plot of a 10^ point distribution on the unit square generated by Model 6. 
G. 2-D Unstructured (Model 7) 

Looking at Fig. □ (Model 2), Fig. [TT] (Model 4) and Fig. [H (Model 6), we notice square 
shaped artifacts and diagonal point formations that arise due to the structured, square-based 
nature of those point placement schemes. Model 7 is a two dimensional point placement 
scheme that aims to avoid such artifacts, while maintaining the two ingredients of reproduc- 
tion and local resettlement. We begin with a single point and assign to this point a child 
distance dc = 1/4:. At each subsequent discrete time step t, each preexisting point (parent) 
is selected to "reproduce" with probability p, where the reproduction consists of a child 
point being placed at a random location on the circumference of a circle centered on the 
parent point and with radius equal to the parent point's dc- A child distance equal to half 
the parent point's dc is assigned to the child point. 
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Model 7 results in multifractal point distributions for all values of p. This is in contrast to 
the trivial dimension spectrum Dq = 1 for the p = 1 limit of Model 1 and the trivial Dq = 2 
result for the p = 1 limit of Models 2 and 3. In addition, while all the models presented 
thus far possess an integer valued box counting dimension, Dq, equal to the dimensionality 
of their respective ambient spaces, Dq ^ 2 for all p in Model 7. 

In a manner similar to the p ^ limits of Models 1 and 2, we can define the p ^ 
limit of Model 7 and create our point distribution one point at a time. The construction 
is as follows: At each time step a parent point is chosen from all preexisting points with 
equal probability. The parent point then produces a child according to the reproductive 
step described above in the general p case. A distribution generated by the p — > limit 
of Model 7 is shown in Fig. ^| As intended, there are no grid-type artifacts visible. The 
p —>■ limit of Model 7 produces point distributions with Di = 1.36 ± 0.05 (averaged over 
20 distributions of 10^ points). 

H. Discussion 

Table m summarizes our results for Di for our models and for the Earth at Night. Ta- 
ble In] summarizes gives qualitative descriptions of the reproduction and local resettlement 
mechanisms of our seven models. 

Our result that the different strategies of reproduction and local resettlement used in 
these models f Table ITT]) all lead to multifractal spatial patterns of population density strongly 
suggests that multifractality may be a generic feature in real situations in which processes 
involving reproduction and local resettlement take place. 

We emphasize that the spirit of our approach is very minimalist. In many real cases 
multiple interacting complex processes undoubtedly influence the determination of popu- 
lation patterns. For example, several important factors for the EaN image are geography 
(mountains, rivers, deserts, etc.), politics, societal and cultural factors, economics, etc. Our 
models show that considerations of underlying spatial heterogeneities are not required for 
explaining the existence of fractally heterogeneous distributions: even very simple dynamics 
incorporating reproduction and local resettlement are sufficient. 
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FIG. 15: A plot of a 10^ point distribution on the unit square generated by the p ^ construction 
of Model 7. 

III. THE EFFECT OF CLIPPED DATA ON MULTIFRACTALITY 

We now revisit the North American part of the Earth at Night image (Fig. ^ and analyze 
its fractal properties. Figure El is a plot of Dq calculated from the measured light distri- 
bution, where the error bars indicate the uncertainty involved in extracting the slopes of 
the various scaling regions. Although we observe a multifractal Dq curve, studies on human 
population distributions 0, [l^ suggest that the a steeper Dq should result. We now show 
how this discrepancy can be resolved. 

Figure El is a log- log plot of the histogram of pixel intensities for North America in the 
EaN image. A spike is observed for a range of intensity values near the maximum value. 
This suggests that some regions with high light intensities caused some of the photoelectric 
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TABLE I: A summary of numerical and theoretical results for the point placement models and 



the Earth at Night image (Fig. . 
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FIG. 16: Plot of Dq for North America in the Earth at Night image. The error bars reflect the 
uncertainty involved in determining the slopes of the scaling regions in the plots of the quantity 
(1 — q)^^ ln/((7, e) versus ln(l/e) [see Eq. 
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TABLE II: Summary of reproduction and local resettlement strategies employed in Models 1-7. 



Model 


Reproduction 


Resettlement 


1 


parent has one 
child 


child settles in the middle of 
chosen interval 


2 


parent reproduces 
by four-fold 
cell division 


children share equally the 
original square area occupied by 
their parent 


3 


parent has 
three children 


parent remains in place; 
children settle equal shares 
of parent's triangular area 




or 




parent reproduces 
by four-fold 
cell division 


children share equally the 

original triangular area occupied by 

their parent 


4 


parent has one 
child 


child settles randomly chosen 
adjacent square 


5 


parent has one 
child 


child settles largest 
adjacent interval 


6 


parent has one 
child 


child settles largest and 
sparsest adjacent square 


7 


parent has one 
child 


child settles 

in a random direction 



cells of the satellites' sensors to saturate, thereby clipping the intensity values at a maximum 
allowable value. Additionally, it is known that, when a photoelectric cell is subjected to high 
intensities, it can trigger surrounding cells to register more light, resulting in the so called 
"blooming" effect, which could account for the broadness of the observed high-intensity 
spike. We note that the saturated pixels make up about 8% of the total number of non-zero 
intensity pixels. We expect that blooming only changes intensities by a marginal amount, 
whereas chpping can make a drastic change. Thus the effect of blooming on Dq should be 
negligible and we ignore it. 
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FIG. 17: Histogram of individual pixel intensities for North America in the Earth at Night image. 

To illustrate the effects that this type of clipping has on multifractal distributions, we 
apply a scheme that mimics the chpping in the EaN image to a point distribution generated 
by one of our models. For this purpose we choose a 3 x 10^ point distribution generated by 
the p ^ limit of Model 2. (Similar results are obtained using our other models.) First we 
roughly match the smallest scale of our distribution to that of the EaN image. Since the part 
of the EaN image we analyze is made up of 3000 x 6000 pixels, we cover our distribution with 
a fixed size grid of simulated pixels in which a simulated pixel corresponds to a grid square 
with edge length of 2~^^, which is about 1/V3000 x 6000. Thus, while each pixel of the EaN 
image represents a light intensity value, each simulated pixel contains a number of individual 
points. The clipping procedure involves picking a maximum number of points allowed in a 
simulated pixel and changing the number of points in any simulated pixel containing more 
than this amount to the maximum amount. In order for the clipped point distribution to 
exhibit an appreciable range for the numbers of points in the individual simulated pixels 
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a) Number of points in a pixel b) Number of points in a pixel 



FIG. 18: Histograms of the number of points in individual simulated pixels (boxes in an e = 1/4096 
grid) used to cover a 3 x 10^ point distribution generated by the p — > limit of Model 2, a) for the 
original distribution and b) after the original distribution was clipped by allowing a maximum of 
40 points in any simulated pixel. 

we choose the maximum clipping value such that the total number of clipped simulated 
pixels make up only 4% of the total nonempty simulated pixels (this is in contrast to the 8% 
observed in the EaN image). For this particular distribution, any simulated pixel containing 
more than 40 points has its count reset to 40. 

Figures ITHT a) and llSf b) are histograms showing the distribution of the number of points 
contained in individual simulated pixels before and after the clipping is applied, respectively. 
We see that the clipped histogram is qualitatively similar to the histogram for the EaN image 
(Fig. ITTj) . Next, we calculate Dq for the point distribution before and after clipping (open 
circles and open triangles in Fig. [THl respectively). We see that the clipping procedure 
transforms the undipped Dq to a flatter curve, making it more similar to the Dq for North 
America in the EaN image (Fig. IT^ . Thus, we are able to conclude that the EaN image 
characteristics are consistent with a multifractal light intensity distribution with a steep Dq 
sensed by an instrument that saturates at a maximum measurable intensity value. 
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FIG. 19: Plot of Dq for a 3 x 10^ point distribution generated by the p — > limit of Model 2 before 
the clipping procedure discussed in Section ITTTI is applied (open circles), and after it is applied 
(open triangles). The error bars reflect the uncertainty involved in determining the slopes of the 
scaling regions in the plots of the quantity (1 — q)^^ ln/((7, e) versus ln(l/e) [see Eq. (0)]. 

IV. INHOMOGENEITY 

Our models can be generalized to include geographical inhomogeneity. For example, 
one can imagine that the underlying space is supplemented by a fertility field, and that 
reproduction is more likely in regions of higher fertility. To investigate the effects of such 
inhomogeneities on the Dg of a point distribution, we modify the p —>■ construction of 
Model 2, by use of a fertility field F{x,y) which we use to construct a space-dependent, 
parent node selection probability. If the location of point i is {xi,yi), then the selection 
probability for this point is taken to be Pi = F{xi,yi)/ J2j F{xj,yj) (this is in contrast to 
the original prescription of choosing the parent point from all preexisting points with equal 
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probability). We investigate two forms for the fertility field. First we use a smooth field 

Fs{x,y) = 1 - 0.3 {cos(27rx) +sin(27ry)} (4) 

on the unit square (0 < a; < 1, < ?/ < 1). We generate point distributions with this 
modified version of the p —>■ limit of Model 2 and find that, while the regions of larger 
Fg in the center of the square are much more dense with points than the regions of smaller 
Fg near the corners, there is, nevertheless, no discernible difference in Dg compared to the 
homogeneous case. That is, both the values of Dg, as well as the quality and extent of the 
scaling ranges, remain unchanged. 

Next, to explore whether the above result can be attributed to the local smoothness of 
Fg, we consider a rough fertility field, 

Fr{x,y) = l + 0.1[Ux) + Uy)], (5) 

oo 

/^(z) = -^a-^cos(27r/3^z). (6) 

j=0 

with a = 1.5 and P = 3. For a < P, the function fw{z) is a 'Weierstrass function' (see 
Fig. 120)); it is rough in the sense that, although it is continuous, it is nondifferentiable, 
and the graph of fw{z) versus 2; is a fractal curve (the fractal dimension of the curve is 
2 — {In a) / {In l3) 1.63 ^|). With a = 1.5, Fr has the same range of variation as Fg 
(i.e., 0.4 to 1.6). Applying this fertility field to the p — limit of Model 2 results in the 
same values and scaling ranges for Dg as in the homogeneous case and the smooth fertility 
field case [Eq. (0}]. This suggests that, while inhomogeneities of the underlying space may 
dictate certain aspects of the distribution of a growing population, for example, where the 
points are more likely to settle, its fractality can be due mainly to reproduction and local 
resettlement processes. 



V. THEORY 
A. Model 1 



Here we derive an analytical expression for the fractal dimension spectrum Dg of popu- 
lation distributions generated by the one dimensional point placement scheme of Model 1. 

We define N{k, t) to be the average (over different realizations of Model 1) number of 
intervals of length at time t. The initial condition for this model (at t = 1) is a single 
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FIG. 20: A plot of the Weierstrass function fwi^) [Eq. In our numerical implementation the 
infinite upper limit on the summation in Eq. Q is replaced by jmax, where jmax is chosen such 
that the smallest value of e used in determining Dq is at least one order of magnitude larger than 
the smallest scale of the roughness, 3~-?™-"^ . 

point on the unit length circumference of a circle, implying N{k, 1) = 6ko, where 6kj is the 
Kronecker delta. At each time t, each preexisting interval is bisected with probability p. 
Thus, with probability p an interval of length is replaced by two intervals of lengths 
2^('^'+i). The discrete time evolution equation for is then 



where pN{k — 1, t) is the average (over different realizations of Model 1) number of intervals 
of length 2^*^'^'^^'' bisected at time t + 1, and (1 — p)N(k,t) is the average number of 2"*^ 
length intervals that remain unchanged at t + 1. We now define N{t) = 1]^ A^(A;,t) to be 
the average number of intervals (of all sizes) at time t. Taking the sum over k of Eq. ^ we 
obtain the simple recursion relation 



which, since A^(l) = 1, gives the result N{t) = (1 + for the evolution of the average 
total number of intervals. 

In Appendix El we introduce the partition function formalism for calculating Dq. We 



N{k, t+l) = 2pN{k - 1, t) + (1 - p)N{k, t) 



(7) 



N{t + l) = {l+p)N{t), 



(8) 
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define the quantity t = {q — l)Dq and derive the following expression for g(r) for Model 1 
[Eq. (IXTO)) ]: 



g(r)-l + ^hm j-^, (2 ), - ^^^^^ . (9) 

Multiplying both sides of Eq. (0) by 2^'^ and then summing over k, we obtain the recursion 
relation, 

= {2'2p + 1 - p) {2'^)„ (10) 



(2^0t= fT^(2^2j9+l-J9)) . (11) 



which, since (2 = 1, gives, 

r2^2n + 1 - nH 

1 +p 

Substituting Eq. (fTT|) into Eq. Q, we obtain 

In (2^2p + l-p) 



(12) 



ln(l+p) 

Then, since Dg = T{q)/{q — 1), we invert Eq. (fT^ for r(g) and obtain the following result 



for Dg-. 



_ log, (£^1± (1±£) 

In particular, using L'Hospital's Rule 



= ^ (13) 



^ ^ (l+p)log.(l+p) . 

^ 2p ■ ^ ^ 

For p = 1 we get the expected result Dg = 1, while for the ]? — > limit, we obtain 

log2(l + g) - 1 , . 

Dg = ^-^ (15) 

and 

fl, = ^.0.72. (16) 



B. Models 2 and 3 



Here we derive the analytical Dg for Models 2 and 3. We adopt the area partitioning 
viewpoint presented in Sec. IIIBl where we focus on the squares containing the points rather 
than the points themselves. We begin by defining N{k, t) to be the average number of 
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squares with edge length at time t. The initial state (at t = 1) is a single square (with 
a point at its center) with unit edge length, giving us N{k, 1) = 6ko. At each time t, each 
preexisting square is divided into four equal sized squares with half the edge length of the 
original square with probability p (Fig. E)). Thus, with probability p a square of edge length 
is replaced by four squares of edge lengths 2"*^*^"'"^''. The discrete time evolution equation 
for N{k, t) is then 



We define N{t) = 1]^ A^(A;, t) to be the average number of squares (of all sizes) at time t, 
and take the sum over k of Eq. (fT7|) to obtain 



which, since A^(l) = 1, gives the result N{t) = (1 + SpY^^ for the evolution of the average 
total number of intervals. Interchanging the squares in the above derivation with equilateral 
triangles, one can see that the same evolution equation holds for the triangles in Model 3. 

In Appendix ^ we show that Eq. ©, which allowed us to calculate g(r) in terms of 
{2^'^)t for Model 1, holds for Models 2 and 3 as well. Following the procedure in Sec. IV A| 
we multiply both sides of Eq. (fTTj) by 2^^"^ and sum over k, solving the resulting recursion 
relation to obtain 



N{k, t + l) = ApN{k - 1, t) + (1 - p)N{k, t). 



(17) 



N{t + l) = {l + 3p)N{t), 



(18) 




(19) 



Substituting Eq. (|TT?|) into Eq. ©, we obtain 

_ In (2^4p +l-p) 
^ ~ ln(l + 3p) 

Thus, for Models 2 and 3 we obtain the multifractal dimension spectrum 



(20) 




(21) 



- 1 



and the information dimension 




(23) 



(22) 



and 



3 



(24) 



4 log 2 



1.08. 
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C. Model 4 



The calculation of Dq for Model 4 is based on the classification of squares into types I 
and II, discussed in Sec. IllFI We begin by defining Ni{k,t) and N2{k,t) to be the average 
(over different realizations) number of type I and type II squares, respectively, with edge 
length 2~'^ at time t. The initial state (t = 1) is a single point at (0.5,0.5). The point 
divides the unit square into four equal squares of edge length 1/2, with each square having 
one point on one of its vertices. Thus Ni{k,l) = 45^1, while N2{k,l) = for all k. At 
each time t, the new point is equally likely to appear in each type I square, but is twice 
as likely to appear in each type II square, because there are two parent points that can 
produce offspring in each type II square. In other words, a type I square with edge length 
2^'^ is chosen with probability Ni{k,t)/ J2j i^iU^t) + 2A^2(j, i)), while a type II square of 
the same size is chosen with probability 2N2{k,t)/ J2j {Ni{j,t) + 2A^2(j, ^))- When a type I 
square is chosen, the new point is placed in its center, destroying the original square and 
creating three type I squares and one type II square, all with half the edge length of the 
original square. On the other hand, when the new point is in a type II square, it is replaced 
by two type I squares and two type II squares, all having half the edge length of the original 
square. Notice that either way, the quantity J2j i^iU^ t) + ^^^2^, t)) increases by 4 between 
time t and t + 1, so that J2j t) + 2N2{i, t)) = At. 

The discrete time evolution equations for both types of squares are then 

iVi(fc, t + l) = N,{k, t) (l - ^) + Ni{k - 1, t)^ + mik - 1, t)^, (25) 

N2ik, t + l)= N2ik, t) (l - ^) + Nr{k - 1, t)^ + 2iV2(A; - 1, t)^. (26) 
In Appendix El we show that for Model 4 [Eq. ()A14|l ]. 

ln(Ti(t)+2^T2(t)) 
= h^^ ' (27) 

where Ti{t) = Y.k'^^^ Ni{k,t). Multiplying Eqs. ^ and ^ by 2^^ and summing over /c, 
we obtain 

m + 1) = T,{t) (l + ^-^) + Ut)^-^ (28) 

T2{t + 1) = T2{t) (l + ^^^^) + Ti(t)|. (29) 
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At this point, we make the continuous time approximation Ti{t + 1) — Ti{t) ^ dTi{t)/dt, 
vahd for large t, which allows us to represent Eqs. ()28|) and as the matrix differential 
equation, 



where 



dT(.) 
ds 



T{s) 



M-T(s), 



T2is) 



(30) 



(31) 



M 



(32) 



3-2^-1 4-2^ 
2^ 4-2^-2 
and s = (1/4) Int. 

Equation pn|l is solved in terms of the eigenvalues Xa, Xb and the eigenvectors T^, T;, of 
the matrix M, giving 



(33) 



for i = a,b, with A^, A;, distinct and real and A^ > Xb- We can express the quantity (Ti(t) + 
2'^T2{t)) in Eq. in terms of a linear combination of the components of these eigenvectors, 
whose long time behavior is dominated by the behavior of T^, the eigenvector associated 
with the larger eigenvalue A^. Thus (|77j) gives us the simple result 



q = A./4. 

The eigenvalues of M are the roots of its characteristic polynomial. 



A^ + A(3 - 7 ■ 2^) + 2 - 10 • 2^ + 8 • 2 



2r 



0. 



This equation is also quadratic in 2"^, allowing us to solve for r(A) 



r(A) = log2 



— (10 + 7A ± \/l7A2 + 44A + 36 
16 ^ 



(34) 



(35) 



(36) 



We substitute the relationship from Eq. into Eq. and, since r(g) = {q — 1)-Dq, we 
resolve the ambiguity of the ± sign by requiring that r(g = 1) = 0. This results in 



Dq = logs 

q-1 



5 + 14g - A/68g2 + 44^ + 9 



(37) 



for the dimension spectrum of Model 4. In particular, applying L'Hospital's Rule to (|T7j). 
the information dimension is 



11 In 2 



1.05. 



(38) 
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VI. CONCLUSION 



The main conclusion of this paper is that reproduction and local resettlement processes 
may lead to multifractal spatial distributions for growing populations. We introduced a 
number of point placement models in one and two dimensions and showed that the models 
resulted in multifractal distributions. Furthermore, we have demonstrated qualitative simi- 
larity between the example of the Earth at Night image and clipped versions of distributions 
generated by our models. We thus suggest that the mechanism by which our models create 
multifractal distributions may be operative in the growth of real systems. 

We thank Dan Lathrop and Rowena Ball for useful conversations. This work was sup- 
ported by ONR (Physics) and by NSF (Contract Nos. PHYS 0098632 and DMS 0104087). 



APPENDIX A: CALCULATION OF MODEL DIMENSION SPECTRA 
1. The Partition Function Formalism 

The partition function formalism ^ |^ is an alternative to the fixed sized grid method 
of Eq. for calculating a dimension spectrum of a measure. We demonstrate here how 
it allows us to relate the various interval and square size distribution functions found in 
Models 1, 2, 3 and 4 to the dimensions of the point distributions generated by them. 

We cover the measure of interest with a disjoint covering {Si}, i = 1,2...N, where each 
element Si of the covering set has a diameter less than or equal to 6. (The diameter is 
the largest possible distance between two points in 5*^.) The partition function is defined as 



N 



r,(r,{5,},5) = E/ifAI, (Al) 

4 = 1 

where /ij is the measure of Si. For a given 6 the covering {Si} is now chosen such that 
Eq. (jAljl is maximized (for g > 1) or minimized (for g < 1), which defines 

mfs^Tg{T,{Si},5), for g < 1 . 



Then, letting 5 — > 0, we define 



r,(r) = limr,(r,<5). (A3) 

0— +U 
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The quantity Tq{T) experiences a jump from to +00, as r is increased, at a critical value 
that we denote T(g). The dimension Dq of the measure is then defined as 

D, = Tiq)/{q-l). (A4) 

In practice, it is difficult to determine whether a particular covering {Si} is optimal in 
the sense of Eq. (IA2|) . However, one can often compute the correct dimension by consid- 
ering specific coverings, in the following sense. Consider a sequence of disjoint coverings 
{S^^^} where covering m of the sequence has maximum diameter (eS""^ < for every 
component i of covering m), and S^"^^ converges to zero as m approaches infinity. If the limit 

f,(r)= hm r,(r,{5^)},5(™)) (A5) 

exists, then like Tq{T) it experiences a jump from to +00, as r increases, at a critical value 
r = 'r(g). In terms of the sequence {S^"^^} of coverings (not necessarily optimal), we can 
then compute a dimension spectrum 

Dq = fiq)/iq~l), (A6) 

Since the sequence of coverings used to compute Dq may be suboptimal, we have f(g) > 
r(g) for g > 1, and f{q) < T[q) for g < 1. Thus Dq > Dq in each case. Assuming the 
limit in Eq. that defines Dq to exist, then Dq can be computed by the partition function 
formalism above using equal size coverings 

tM) = for every i and m). Thus Dq > Dq 
as well. While there do exist examples for which Dq^ Dqioi q <1 \^h\ the two dimensions 
are typically found to coincide in analytical examples with physical bases. Furthermore, for 
q> 1 the two dimensions always coincide: 

Theorem: For every probability measure /i for which the limit in Eq. exists, Dq = Dq 
for all g > 1. (If one defines Dq as a liminf, then this equality holds for all /i.) 
We prove this theorem in Appendix IbI 

The theorem demonstrates that, for g > 1 at least, the dimension spectrum computed 
according to the partition function formalism from a particular sequence of coverings is not 
that sensitive to the type of covering; equal size coverings yield the same spectrum as a 
sequence of coverings that is optimal in the sense of Eq. ()A2|) . While it may be possible 
to compute a different value of Dq from a sufficiently suboptimal sequence of coverings, we 
conjecture that for the coverings we consider below, 

Dq = Dq = Dq. (A7) 
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Furthermore, we conjecture that Eq. ()A7|) holds not only for g > 1, but for all g > 0. 

For our purposes, the application of Eq. ()A7|) for g < 1 is supported by the following 
reasoning. In all of the cases below, the formula we obtain for Dq is an analytic function 
of q (for our region of interest < g < 2). If Dg = Dg for g > 1, it follows that if Dg is 
continuous at g = 1, then Di = Di. Furthermore, if Dg is analytic, then Dg = Dg for g < 1 
as well. 



2. Application to the Models 

The practical implication of Eq. ()A7|) for us is significant. It means that we can choose 
any sequence of coverings whose maximum diameter converges to zero, and this choice can 
be made in a manner that facilitates analytic computations. The result will be the same as 
that for an equal cube size covering [Eq. (^] or an optimal covering [Eq. ()A2|) ]. 

For Model 1, we choose the intervals between the points at time t as a covering for the 
distribution. We regard each interval as covering the equivalent of one point since it contains 
two half points at each of its ends. Hence, on average, /Xj = {J2k^{k,t))^^ = 1/N(t) for all 
i, and Eq. ()A5|) becomes [l^ 

f,(r) = lim(l/Ar(t)r^er, (A8) 

i 

where the 5 — limit is replaced by t oo, since, as t increases the size of the largest 
interval decreases to zero with probability one. Since Inf^(r) = — oo for r < f{q) and 
Inrq(r) = oo for r > r(g), we equate InTglr) to zero and obtain g as a function of the 
transition value f, 

,(^, . = '"g'^'f-f". (A9) 

We then have [see Eq. Q] 

where 

Each square in Model 2 and triangle in Model 3 contains, on average, {J2k ^{k,t))^^ = 
1/N{t) of the total measure. Thus, like we derived for Model 1, we obtain 

g(r) = 1 + lim ^^^^TT- (A12) 
^ t^^ IniV(t) ^ ^ 
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For Model 4 we use the two types of squares (type I and II) to cover the distribution. A 
type I square has one point on one of its vertices and hence covers the equivalent of a quarter 
of a point, while a type II square, with two points on its vertices, covers the equivalent of 
half a point. Thus, since at time t we have t total points, the measure contained in a type I 
square is l/4t, while that contained in a type II square is l/2t. Thus Eq. ()A5|) becomes 

r'M^.'^iio^g'r + a'g^rJ, (A13) 

where the first summation is over all type I squares and the second is over all type II squares. 
Taking the logarithm of this expression and equating it to zero, we obtain 

where Ti{t) =j:f,2''^Ni{k,t). 



APPENDIX B: PROOF OF THEOREM FROM APPENDIX El 

As we noted in Appendix \^ Dq < Dg for all a > 0, so we now show that also Dg > Dq 
for g > 1. Our proof is based on the theorem in [l^ that for g > 1, the quantity 




is finite for s < Dq and infinite for s > Dq. To see that Dq > Dq, we show that Dq > s 
whenever Iq{s) is finite. Let {Si\ be a disjoint covering of as in Appendix |Xl and recall 
that /ij is the measure of Si and is its diameter. Then considering only the contribution 
to Iq{s) from points x and y that are in the same Si, we have 




N 
i=l 

= Tq{{q-l)s,{S,},6) 
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by ()A1|) . This implies that Tq {{q — l)s,S) < Iq{s) by ()A2|1 . since g > 1, and hence 
Tq ((g — l)s) < Iq{s) by ()A3|) . Therefore if is finite, then ((g — l)s) is finite, whence 
r{q) > {q- and finally Dq > s by (IX4|l . 
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where A(e) is the number of e-intervals needed to cover the set. Next we observe that the 
distance between 1/n and l/(n + 1) is approximately 1/n^ for large n. Thus setting e = 1/n^, 
we need n intervals to cover the first n elements of the set, one interval for each such element. 
To cover the remaining elements we must cover the interval {l/{n + 1),0). This requires 
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(l/e)[l/(n + 1)] = n^/{n + 1) ~ n intervals. Thus iV(e) ~ 2n = 2/6^2, yielding Dq = 1/2. 
[16] There is a subtlety here. As t increases we are changing both the covering and the distribution of 
points. We really should be looking at the limiting distribution and not the time t distribution. 
Since each interval is equally likely to subdivide, the expected measure in each of the time t 
intervals is the same, and our computations in this section reflect this approximation. 
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